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Abstract 

An ultracold Fermi atomic gas at unitarity presents universal properties that in the dilute limit 
can be well described by a contact interaction. By employing a guiding function with correct 
boundary conditions and making simple modifications to the sampling procedure we are able to 
calculate the properties of a true contact interaction with the diffusion Monte Carlo method. The 
results are obtained with small variances. Our calculations for the Bertsch and contact parameters 
are in excellent agreement with published experiments. The possibility of using a more faithful 
description of ultracold atomic gases can help uncover additional features of ultracold atomic 
gases. In addition, this work paves the way to perform quantum Monte Carlo calculations for 
other systems interacting with contact interactions, where the description using potentials with 
finite effective range might not be accurate. 

PACS numbers: 67.85.Be, 03.75.Ss 
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I. INTRODUCTION 


Systems formed by fermions have many-body properties that are of central importance for 
understanding observed phenomena in many helds of physics. These helds include ultracold 
gases, condensed matter, and nuclear physics. The possibility of handling ultracold atomic 
Fermi gases, in a very precise way, has allowed testing quantum many-body theories in an 
unprecedented set of experimental conditions. 

Ultracold Fermi gases can be tuned from weakly interacting to strongly correlated regime 
by application of magnetic helds near a Feshbach resonance [1]. When the interaction has 
diverging scattering length, the unitary limit, the system presents universal properties, i.e., 
it does not depend on the nature of the interactions. The system universality allows one 
to study the crossover from the Bardeen-Cooper-Schrieffer (BCS) superhuid state to the 
Bose-Einstein condensed (BEC) state, in general [2]. 

Countless efforts were made and continue to be made [3] to uncover the many aspects 
involved in the observed phenomena presented by the ultracold Fermi gases. In the present 
work we investigate the unitary limit of this system at the crossover from BCS to the BEC 
regime with an s-wave contact interaction. 

Interactions of two neutral atoms are not always easy to describe accurately. However in 
the dilute regime, interactions can be well represented by two-body collisions using a con¬ 
tact potential. Nevertheless a straightforward consideration of this type of potential makes 
theoretical investigations problematic because when two particles approach one another the 
wave function diverges. This difficulty is usually avoided by adopting pseudopotentials of 
the Poschl-Teller, hard sphere, square well, or other forms [4]. In this fashion, valuable 
insights have come from quantum Monte Carlo methods [5-8], despite the fact that hnite- 
range potentials lead to incorrect scattering properties, which are fundamental quantities of 
these systems. The resulting calculations must therefore include an additional extrapolation 
to zero range. Since the trial wave functions diverge in this limit, the extrapolations are not 
well behaved in this limit. 

It is not just a matter of principle or of interest in itself to avoid using hnite-range 
pseudopotentials to describe the two-body interaction of ultracold Fermi gases. For instance, 
it is important to avoid the possible influence of the true ground state of the Poschl-Teller 
model system, since it may have tightly-bound states highly dependent on the chosen range. 
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For repulsive interactions, there are still open questions about the ferromagnetic character of 
the ground state and what kind of ferromagnetic transition the system undergoes in this case 
[9-11]. The possibility of simulating Fermi atomic gases considering a contact interaction 
will help solve questions like those previously mentioned. On the other hand, studies of Bose 
systems, including Bose-Fermi mixtures, have been done only using hnite range interactions 
in quantum Monte Carlo calculations, see for example Refs. [12, 13], introducing possible 
bias in the calculation. 

The contact interaction as we have considered allows the quantities of interest to be 
obtained without the additional burden of performing extrapolations to zero-range interac¬ 
tions. This is useful in a twofold way. It can help understand how previous results might 
have been affected by the use of hnite range potentials, and also because the calculations 
become simpler. Moreover, the results we present depend on relatively small changes of the 
standard diffusion Monte Carlo (DMC) algorithm. Additionally, we show how to compute 
the two-body propagator for particles interacting through a contact potential, which is an 
interesting result in itself. 


II. METHODOLOGY 


The system we study consists of N fermionic particles described through the Hamiltonian 

( 1 ) 


H = - 

2m 


N/2 N'/2 




where the hrst terms are the kinetic energies of the up-spin (unprimed index) and down-spin 
(primed index) particles and the last term is the zero-range interatomic potential. Here we 
focus on the unpolarized system, and N/2 particles are spin-up, and N/2 are spin-down. The 
simplest way the solve Eq. (1) is to introduce a trial variational wave function where 

R = {ri, r'^, ■ ■ ■ , rp^/ 2 , r 7 V 72 }, and minimizing the expectation value of H [14]. Typically, one 
samples M conhgurations from the probability density proportional to and average 

the local energy. The value of the variational Monte Carlo energy Evmc is a ground state 
upper bound and it normally depends on the quality of the trial wave function. Beyond 
the VMC calculation using the diffusion Monte Carlo(DMC) one can project out the lowest 
state of the system from 
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The Schrodinger equation can be written as a diffusion equation in imaginary time r = 

it/h, 

( 2 ) 

where Et is introduced to stabilize the norm of ip{R-,T) in the limit of r —)■ cx). It is 
convenient to expand the trial wave function in the eigenstates {ipi} of H and make Et ~ Eq. 
In that case, it is straightforward to show [15, 16] that we project out the 'ipo evolving 
in the imaginary time. In practice, the total Monte Carlo time is divided into n small time 
steps such that r = At ■ n and, the exact wave function is propagated like 

ij{R-,T)= lim '^t{R) . (3) 

n^oo *- 

In the DMC method the projected state at large imaginary time is the lowest energy state 
not orthogonal to the trial function. 

The zero-range interaction in s-wave is enforced by using an importance function that 
satishes the Bethe-Peierls condition (l/rjj/ — 1/a) when —)■ 0, where a is the two-body 

scattering length. At unitarity, (a —)■ cx)) the Bethe-Peierls condition reduces simply to the 
wave funtion being proportional to the inverse modulus of the pair relative distance at small 
separations. This approach allows us to treat the system with the zero-range pseudopotential 
as formed from pairs of free particles subject to the correct boundary conditions when a pair 
separation distance goes to zero. 

The importance function we use can be written as: 

'J'r = ; (4) 

ij' 

^Bcs = -4, [0(rii/)0(r22') ■ ■ • 0(rAr/2 Af 72 )] (5) 

where the Jastrow pair function /(r) correlates the unlike spin pairs of the system. The pair 
function is chosen to satisfy the Bethe-Peierls condition, and the antisymmetric BCS function 
^Bcs is well behaved at small pair separations and dehnes the nodal surface structure [5]. 
The operator A antisymmetrizes like spin pairs. Here we use the general form for the BCS 
part where the pairing functions are written like a set of plane waves 

Ns 

= ( 6 ) 

i=i 

and Cj are variational parameters. The Cj coefficients with the same magnitude kj are equal. 
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These orbitals have the same form of previous works [16]. The short-range pairing func¬ 
tion of Refs. 5 and 8 is not included since the boundary condition from the potential is 
enforced by the Jastrow factor. We have considered Ns = 20 shells in the guiding function 
and have obtained converged energies for both variational and DMC calculations. The co¬ 
efficients entering in the pairing orbitals have been optimized as described in [8] using the 
stochastic reconhguration method [17]. The function ^bcs is projected on the subspace 
with hxed number of particles N 


\BCS) = -f 


( 7 ) 


If for ki > kp all the n* are equal to zero this function reduces to a product of two Slater 
determinants of plane waves[16]. We call this the Jastrow-Slater wave function; it was used 
in our previous work [16]. 

The Jastrow pair function /(r) in Eq. (4) correlates the unlike spin pairs, and we take 

d cosh(Ar) 


f{r) = 


( 8 ) 


r cosh(A(i) 

with f{r>d) = l, the parameter A is chosen to make / and its hrst derivative continuous at 
the healing distance r = d. Its value is of the order of the inverse interparticle distance and 
is determined with a variational calculation. It is important that the Jastrow pair function 
has the correct boundary condition at short distances. 

Our variational calculations are performed as in Ref. 16. The variance of the energy of 
this trial function with the usual VMC method is not well behaved. This can be seen by 
looking at the form of the trial function when a pair separation is small. For example, if 
up-spin particle 1 and down-spin particle 1' are close together, the Jastrow factor /(rn/) 
goes like The term in the Hamiltonian where — -|- V^/) operates on this gives 

zero except at the origin where it cancels the contact interaction, so that part of the local 
energy is well behaved. The problem terms are those like Vi/(rii/) ■ Vi This term 

goes like ■ Vi at small distances. The r?,, term in the volume element as well as 
the angular integration makes this term give a well behaved contribution to the energy as 
the pair separation goes to zero, however, this term is squared in the energy variance, and 
the variance diverges. 

The exact ground state wave function must, of course, have zero variance with = 

Eq independent of R. This means that the exact wave function must have additional terms 
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which cancel these divergences. These would have the form of three-body correlations which 
would cancel the divergences from terms like Vi/(rii/) ■ Vi/(rij/), and backflow terms to 
cancel divergences from terms like Vi/(rii/) ■ Vi(;/)(rij/). Such terms would have to be 
constructed in such a way as to not spoil the necessary boundary conditions. A hierarchy 
of such terms may be required to obtain a well behaved variance of the local energy. 

Since the integrated energy is well behaved, we have chosen to attack the problem by 
modifying the sampling to control the variance. The key insight is that for rn/ —)■ 0, 
interchanging the positions of particles 1 and 1' reverses the sign of the gradient Vi/(rii/), 
but does not change the rest of the trial wave function. We therefore modify the standard 
Metropolis algorithm to include moves which interchange the positions of the closest pair of 
unlike spin particles. If the pair remains the closest pair after interchange, we accept this 
move with the heat bath probability for interchange 


Pint 




int) 


+ 'HURtnt) 


( 9 ) 


where Rint are the coordinates with the closest pair interchanged. We use the method of 
expected values to evaluate the energy, after such a trial move, so that the energy contribu¬ 
tion is EL{Rint)Pint + — Pint)- In the limit of small pair separations, P^t —t and 

the diverging contributions cancel. 

The diffusion Monte Carlo calculations have the same diverging terms in the local energy, 
and we employ a similar technique to control the variance. 

The propagation equation in imaginary time including as an importance function 


^T(R)iiR; r + At) = y <i^R'^^G(R, R'; AT)^T(R')>f’(R', r). (10) 

Since the zero-range interatomic potential is a delta function, the usual Trotter-Suzuki 
decomposition of the propagator is not adequate. The walkers are instead sampled 

from R!] At) [16] and we essentially have only to deal with the kinetic energy 

term of the Hamiltonian. The short time propagator we use is evaluated using the pair 
product form from the two-body propagators g 


G{R',R) = Go{R', fl)n 

Kj 


^o(r',r'.;ri,rj) ’ 


( 11 ) 


where Go and go are the free particle and the free pair propagators, respectively. Note that 
for pairs with the same spin g = go- For pairs with opposite spin g can further be written 
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as ^frei X ^fcm, the product of the relative times the center of mass propagators of the pair 


9(r;,r',;rj,rj.) = a-elW/; r«')gcm(R;j.; R«'). 


(12) 


where Rjj/ = (r* + Vji)/2 is the center of mass of the pair. In our approach it is necessary 
to write the full propagator as above. 

The centers of mass propagate like free particles [16]. On the other hand, the two-body 
propagator is a Green’s function that can be constructed from the normalized solution of 
the of the scattered s-wave function as employed in other papers [16, 18] 

^rei(r, r; Ar) = ^ ipn{r)e~^^^ipl{r') (13) 


where {ip} is the complete set of eigenstates of the two-body Hamiltonian. Since the inter¬ 
action is only in the s-wave, we can separate into partial waves, and the s-wave contribution 
for scattering length a becomes 


9 sir, r',a) 


1 

dTT^rr 


-Re 



(1 ikci') ik{i — r') 

1 + W ^ 


IPk^At 

e ™ -f bound state, 

(14) 


where for positive a, the bound state contribution should be included. The integrals can be 
done straightforwardly in terms of gaussians and error functions, 

1 


9s{r,r',0) = 


gs[r,r ,-oo) = 




STT^rr' V h?At L 
1 


STT^rr' V K^At L 
gsir,r',a) = gsir,r', -oo] 


t L 

TflTl r m f I /\2 m / /\2 

e > +e 


47rrr' a 


At _ r+r 

- 0 ma^ “ erfc 


h?At r + r' I ma^ A 




2a V h?At 


(15) 


where any bound state contribution needs to be added. Here we are primarily interested in 
the unitary case a = —oo where the relative coordinates propagator is particularly simple [19] 

1 


g,e\ir, r'; Ar) = g^^^ir, r'; Ar) + 


-e 4fi^A- 


■ir+r'y 


( 16 ) 


y K^Ar 

where the first term is a free-particle propagator for the relative distances and the last one 
is the contribution of the s-wave scattering. 

The sampling of the importance sampled propagator is accomplished by 

approximating it, sampling the approximation, and using the ratio of to the 

approximation as a weight. 
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We first construct what we call the independent pair propagator Gip{R',R). We sort the 
unlike spin pair distances, and first select the closest pair. We then eliminate all pairs that 
contain the closest pair’s particles. We repeat this process on the remaining pairs. The 
result is a list of independent pairs. The independent pair propagator is the product of the 
pair propagators Eq. 12 for these independent pairs. From the form of the relative pair 
propagator, we see that if the initial separation is much larger than the propagator 

becomes the free particle propagator. Furthermore for large separations, the divergences 
in the trial function can be neglected. Therefore we introduce a cutoff parameter, so that 
if the separation is larger than this parameter, we sample the pair from the free particle 
propagator. If it is less than this cutoff, we approximate the trial function by the Jastrow 
factor for that pair, and for these separations, we take its asymptotic value, given by the 
Bethe-Peierls condition. We sample the center of mass part of the pair propagator from the 
noninteracting center of mass gaussian, and the relative separation from 

ij' 

The details of this sampling are given in the appendix. This general method can be readily 
extended to scattering lengths away from unitarity. 

The value of the cutoff parameter does not affect our results, and for reasonable values 
has very little effect on the variance. We denote the sampled conhguration Ri = R + AR 
obtained as described above by R) = Gip{R, R') ni<j where the product of the 

Bethe-Peirls condition is only over the pairs that are within the cutoff distance. 

To include importance sampling, we use the antithetic “plus-minus” sampling often used 
in nuclear quantum Monte Carlo calculations[20], for the center of mass variables, and the 
relative coordinates beyond the cutoff. For these, the gaussians have the same probability 
of taking the opposite sign. Therefore, it is equally probable for us to have sampled the 
conhguration R2 = R — (i?i — R). 

A divergence in the local energy at the sampled point Ri can occur exactly as in the vari¬ 
ational calculations. To avoid this, two additional conhguration are also considered. These 
are obtained by interchanging the closest pair from Ri and R 2 . Finally, a new conhguration 
is chosen among the Ri according to 


'^g(r;,r) 




(17) 


By performing this choice, the importance sampled ^^^G{R', R) is recovered by through 



the weight 


W{R') 


( 18 ) 




III. RESULTS AND DISCUSSION 

In the unitary limit, the resonant character of the interactions of a Fermi gas makes the 
system have only two possible energy scales: the chemical potential fi and the Fermi energy 
Ep. Therefore these two quantities must be proportional, jj, = ^Ep. As the temperature 
approaches zero, the reduced chemical potential ulEp saturates to the universal value 
Of course, in this limit, fi converges to the system ground state energy. The value of 
^ has been accurately measured: ^ = 0.376(4) [21]. However a more recent work [22] 
suggests corrections to this value, resulting in .^ = 0.370(9). If the atomic interaction 
is described by hnite range pseudopotentials, determining accurate values of ^ requires a 
careful extrapolation to zero effective range [4]. Our result for this quantity, also known as 
the Bertsch parameter, is .^ = 0.390(2), obtained by simulating a system with 66 particles. 
It is obtained in a straightforward manner, subject only to the fixed-node approximation 
and hnite size dependence. The determined value is in reasonable agreement with the 
experimental one. We have observed that there is only a small dependence of this quantity 
on the number of particles in our simulations, as also reported in Ref. [23]. The energy we can 
obtain is in agreement with the best hxed-node diffusion Monte Carlo calculations performed 
using hnite ehective range interactions; in Ref. [24] using the auxiliary-held quantum Monte 
Carlo and a exact lattice technique, ^ was determined as 0.372(5). 

The strong interacting Fermi gases described by contact interactions obey a number of 
universal relations characterized by a single parameter dubbed the contact C by Tan [25]. As 
shown by Zhang and Leggett [26] the contact is able to enclose all of the many-body physics. 
The contact density integrated in the whole space gives the contact, which is proportional 
to the number of pairs with opposite spins that are close together. Its value can also be 
computed straightforwardly from our calculations. Before computing its value it is useful 
to extract a related constant ( from the pair distribution function of unlike-spin pairs as a 
function of the distance presented in Fig. 1. At the unitary limit and for small distances [8] 
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This is because the pair distribution function of particles with opposite spins separated by 
small distances satishes in a hrst approximation g{r) oc /^(r). The behavior of g^^ilkpr) at 
small distances conhrms with what we expect from Eq. (19) as we can verify from the inset 
of Fig. 1. If we modify the £t by imposing 6o = 0 we have estimated ( = 0.755(1). This 
value is slightly smaller that the one obtained with a £t where bo is a free parameter. With 
this fit, it also becomes more clear that a perfect agreement between the DMC results and 
the htted black line in the inset of Fig. 1 occurs for small values of kpr. The BCS result is 
shifted to the right of the Jastrow-Slater, most probably due to a large delocalization of the 
particles in the superfluid state. 

The relation between the constant ( and the contact parameter at unitarity is simple, 
C/kp = 2(/57r [8]. However to make the comparison with experimental data easy, we 
report this quantity in terms of the contact per unit volume given by C/Nkp = 3n'^C/kp. 
Our result, C=2.848(1), is slightly below two recent measurements. A Bragg spectroscopy 
experiment [27] gives the value 3.06 ± 0.08 at the temperature T/Tp = 0.08. A measurement 
using radio-frequency spectroscopy gives 2.9 ± 0.3 at {kpa)~^ = —0.08 and T/Tp = 0.18(2), 
a temperature slightly above the transition temperature Tc [28]. Our computed value is closer 
to the experimental values than previous results determined with a hnite range potential [8]. 

The contact C remarkably controls short-distance correlations. On the other hand, the 
momentum distribution na{k) in the spin state a for large enough momenta is given by 
na{k) = C/k^. We have computed the quantity n{k/kp){k/kpY as a function of k/kp, and 
our results are shown in Fig. 2. The contact term is dominant for momentum states larger 
than approximately l.Qkp, as we can see from the hgure. This dominance is expected since 
n{k/kF){k/kpY —)■ Although the estimated values of n{k/kF){k/kpY are noisy for 

large momenta, it is possible to observe a trend towards the value of C, estimated from the 
pair correlation function, and displayed as a solid line. The less than optimal agreement of 
our results with the experimental data might come from various sources. It might be due 
to calculations done at zero temperature while the experiments are, of course, done at hnite 
temperature. Other possibilities might include the asymptotic form we have considered 
for the guide function; it eventually needs to be improved by including more long-range 
correlations. However it is worthwhile mentioning that other DMC calculations [5] would 
also overestimate the values of n{k/kF) at low values of k. 


10 




kpr 


FIG. 1. (Color online) Pair distribution function of unlike-spin pairs as a function of the distance. 
Our quantum Monte Carlo (QMC) results are for a system with 66 particles. The solid line is the 
best fit of bo + b/{kFr)‘^ to the extrapolated points. For completeness we have included results for 
the Jastrow-Slater model for a non superfluid system of particles [16]. The inset presents the same 
distributions multiplied by (kpr)'^ as a function of the distance. 

IV. SUMMARY 

In summary, we have performed for the first time diffusion Monte Carlo calculations of 
a system interacting with a contact interaction. This has allowed us to have a more faith¬ 
ful description of dilute ultracold Fermi gases at unitarity that opens possibilities of more 
accurate and precise calculations of other important quantities associated with this system. 
The application of this approach has allowed us to compute such quantities as the reduced 
chemical potential and Tan’s contact parameter in better agreement with experiment than 
some previous calculations. We have introduced an alternative way of studying ultracold 
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FIG. 2. (Color online) Momentum density distribution n{k/kF) multiplied by (k/kp)^ as a function 
of k/kp. In the inset we show n{k/kp) as a function of k/kp. The solid line shows our estimated 
value of C. The experimental data is for an inhomogeneous gas [29]. 

atoms in the unitary limit which will be of value in the investigation of these systems, and 
in other situations where a description using a hnite effective range interaction might be 
inaccurate. 
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Appendix A: Sampling the unitary propagator 


We will begin by looking at the dominant part of the wave function when an opposite 
spin pair have a small separation. In this case, we can approximate the trial wave function 
ratio as 

^t(R) . r' 


'hT(-R') r 

Since the propagator consists of the free particle gaussian in all channels except the s- 
wave, we separate it into the usual spherical coordinates r, cos 9, and 0. Starting with the 
free particle propagator, we take the axis along the initial value r' (note we use primed 
coordinate for the initial value here for convenience; in the main text the initial coordinates 
are unprimed and the sampled coordinates are primed), and write the importance sampled 


gaussian as 


Normalizing the cos 9 part 

rl 


r' 1 

e 2 <t2 = 




e 20-2 e ct 2 


cos 9 


, rns 0 2(7 . . / rr 

dcos9e^ =—-sinhf — 


J-i rr' Vo- 

given an r and r' value, we can sample the angular part from 


pM = ^ 


(A2) 


(A3) 


P9{cos9) = ^ , Z 


(A4) 


2a2 sinh(^) 

Once we know r, we can sample cos 9 by sampling a uniform random number 0 < .^ < 1 and 

2 

CT r / 2rr' \ 2rr' 

’ k(l-e-^)+e-^ 


cos 9 = 1 ^ -- In 




The importance sampled gaussian is now 


|r- 


e 2 <t 2 = p^{(j))pg{cos9) 


r (27ro-^)^A 
The integral of the r part over is 




Tiar^ 


e 2 ^ 2—6 2 < t 2 


Hdr-L- 

k spina 


e 2<t 2 — e 2 ct2 


= erf 


a 


V2 


(AS) 


(A6) 


(A7) 


We now look at the “extra” piece from the unitary s-wave interaction. It has the impor¬ 
tance sampled form 


■\/27r _ (i 


r dTT^orr' 


e 20-2 


(AS) 
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Here the angular part is isotropic, so we can sample the angles from 


pM = ^ 

( 0 ) / ^ 
Pe [cose) = - 


and the importance sampled function becomes 

r' ( 0 ), . 

e = P,j>[(p>Pe [cosO) 


r Air'^arr' 

The integral of the r part over r^dr is 


2 1 

-2® 

TT ar^ 


dr\l -e 

TT a 


2 <t 2 = 1 — erf 


a 


V2 


(A9) 


(AlO) 


(All) 


When we add the normalizations of Eqs. A7 and All, we get one, since we should get 
and with the ground-state energy Eq = 0 since there is no bound state for the unitary gas. 

This suggests a way to sample the propagator. We first sample the r value, with probabil¬ 
ity that the r value was sampled from the gaussian we sample cos 6 from pelcos 6), otherwise, 
we sample cos 6 uniformly. In either case, we sample (j) uniformly. 

The basic idea below is to sample from a 1-dimensional gaussian centered at r'. This 
corresponds to the first term of Qs- If the resulting r is greater than zero, then it is a 
legal value. With probability given by the ratio of the radial part of Go divided by the 
sampled 1-dimensional gaussian, we take this r as being sampled from Go, and sample cos 6 

_ (r+r'p 

accordingly. If not, the rejected terms have been sampled from the e 2 <t 2 ^ so they are half 
of the s-wave term. If the resulting sampled r is negative, we take its absolute value and it 
is also sampled from the s-wave term. 

Our algorithm is 

• Sample a 0 uniformly on 0 < 0 < 27r or equivalent. 


• Sample a random variate y from a gaussian with mean zero and variance 1. 

• The r sample is r = \r' + ay\. 

• If (r' -|- ay) < 0 sample cos 9 uniformly. The sampling is complete. 

• If (r' -|- ay) > 0 then sample a random variate 0 < .^ < 1 uniformly. 

_2rr^ 

• If .^ < e sample cos 9 uniformly, else sample cos 9 from p0{cos9). 

A graph of the various terms is shown in Eg. 3. 
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r/a 

FIG. 3. (Color online) The sampling of the radial part of the propagator is illustrated. We sample 
from the gaussian centered around r', in this case taken to be a. The probability of this being from 
the free gaussian is shown as the (red) go region. It is the difference between the sampled gaussian 
and the gaussian centered around —r', also shown. The (blue) gs region is half the probability of 
sampling from the s-wave propagator. If we take the absolute value of the sampled value when 
it is negative, the samples are from the (green) gg w/abs region which gives the other half of the 
s-wave propagator. 
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